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The statistical equilibria of the (conservative) dynamics of the Gross-Pitaevskii Equation (GPE) 
with a finite range of spatial Fourier modes are characterized using a new algorithm, based on a 
stochastically forced Ginzburg-Landau equation (SGLE), that directly generates grand canonical 
distributions. The SGLE-generated distributions are validated against finite-temperature GPE- 
thermalized states and exact (low-temperature) results obtained by steepest descent on the (grand 
canonical) partition function. A standard finite-temperature second-order A-transition is exhibited. 

A new mechanism of GPE thermalization through a direct cascade of energy is found using 
initial conditions with mass and energy distributed at large scales. A long transient with partial 
thermalization at small-scales is observed before the system reaches equilibrium. Vortices are shown 
to disappear as a prelude to final thermalization and their annihilation is related to the contraction of 
vortex rings due to mutual friction. Increasing the amount of dispersion at truncation wavenumber is 
shown to slowdown thermalization and vortex annihilation. A bottleneck that produces spontaneous 
effective self truncation with partial thermalization is characterized in the limit of large dispersive 
effects. 

Metastable counter-flow states, with non-zero values of momentum, are generated using the SGLE 
algorithm. Spontaneous nucleation of vortex ring is observed and the corresponding Arrhenius law is 
characterized. Dynamical counter-flow effects on vortex evolution are investigated using two exact 
solutions of the GPE: traveling vortex rings and a motionless crystal-like lattice of vortex lines. 
Longitudinal effects are produced and measured on the crystal lattice. A dilatation of vortex rings 
is obtained for counter-flows larger than their translational velocity. The vortex ring translational 
velocity has a dependence on temperature that is an order of magnitude above that of the crystal 
lattice, an effect that is related to the presence of finite-amplitude Kelvin waves. This anomalous 
vortex ring velocity is quantitatively reproduced by assuming equipartition of energy of the Kelvin 
waves. Orders of magnitude are given for the predicted effects in weakly interacting Bose-Einstein 
condensates and superfluid 4 He. 
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I. INTRODUCTION 

Finite temperature superfiuids are typically described 
as a mixture of two interpenetrating fluids [1]. At low 
temperatures the normal fluid can be neglected and Lan- 
dau's two-fluids model reduces to the Euler equation for 
an ideal fluid that is irrotational except on (singular) vor- 
tex lines around which the circulation of the velocity is 
quantized. At finite temperature, when both normal fluid 
and superfluid vortices are present (e.g. in the counter- 
flow produced by a heat current) their interaction, called 
"mutual friction" , must also be accounted for [2] . 

In the low-temperature regime the Gross-Pitaevskii 
equation (GPE) (also called the Nonlinear Schrodinger 
Equation) is an alternative description of superfiuids and 
Bose-Einstein Condensates (BEC) [3 |. The GPE is a par- 
tial differential equation (PDE) for a complex wave field 
that is related to the superflow's density and velocity 
by Madelung's transformation [4]. The (non singular) 
nodal lines of the complex wave field correspond to the 
quantum vortices that appear naturally in this model 
with the correct amount of velocity circulation. Just as 
the incompressible Euler equation, the GPE dynamics is 
known to produce [SHE] an energy cascade that leads to 



a Kolmogorov regime with an energy spectrum scaling 
as E(k) <~ fc -5 / 3 . This Kolmogorov regime was also ex- 
perimentally observed in low temperature helium [9l 110] . 
In this experimental context, let us remark that so much 
progress has been made that it is now possible to vi- 
sualize superfluid vortices both in the low-temperature 
regime and in the presence of counter flow by follow- 
ing the trajectories of solid hydrogen tracers in helium 

Several different theories of finite-temperature effects 
in BEC have been proposed and, at the moment, there 
is no consensus on the best model [3] . In one approach it 
has been suggested that, beyond its good description of 
the low-temperature regime, the GPE should also be able 
to describe the classical equilibrium aspects of a finite- 
temperature homogeneous system of ultracold gases, pro- 
vided that that a projection (or truncation) on a finite 
number of Fourier modes is performed . Another 
approach to finite temperatures is the Zaremba-Nikuni- 
Griffin (ZNG) theory [H] which couples the GPE with 
a Boltzmann-like equation for the thermal cloud of non- 
condensed particles. The ZNG theory is known to well 
describe the observed finite temperature decay of soli- 
tons [15] . It also predicts vortex motion in agreement 
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with the standard phenomenology [TB]- In the truncated 
GPE model the small-scales modes are in thermal equi- 
librium. They play the role of the Boltzmann sector of 
the ZNG, somewhat like the (fast) thermalized degrees 
of freedom do in a standard molecular dynamics simula- 
tion. The present paper is devoted to the truncated GPE 
approach. 

Classical truncated systems, that are similar to the 
truncated GPE, have a long history in the context of fluid 
mechanics. Indeed, if the (conservative) Euler equation 
is spectrally truncated, by keeping only a finite number 
of spatial Fourier harmonics, it is well known that it ad- 
mits absolute equilibrium solutions with Gaussian statis- 
tics and equipartition of kinetic energy among all Fourier 
modes Q3H2]]. 

Recently, a series of papers focused on the dynam- 
ics of convergence of the truncated Euler equation to- 
ward the absolute equilibrium. It was found that (long- 
lasting) transient are obtained that are able to mimic 
(irreversible) viscous effects because of the presence of a 
"gas" of partially-thermalized high-wavenumber Fourier 
modes that generates (pseudo) dissipative effects [2"TH2"r?] . 

The main goal of the present paper is to obtain and 
study finite temperature dissipative and counter flow ef- 
fects by extending to the Fourier-truncated GPE the dy- 
namical results that were obtained in the framework of 
the truncated Euler equation. We now give a short re- 
view of what is already known about the truncated GPE 
dynamics. 

The Fourier truncated Gross-Pitaevskii equation was 
first introduced in the context of Bose condensation by 
Davis et al. [T5J as a description of the classical modes 
of a finite-temperature partially-condensed homogeneous 
Bose gas. They considered random initial data defined 
in Fourier space by modes with constant modulus and 
random phases up to some maximum wavenumber (de- 
termined by the energy). They found that, the numeri- 
cal evolution of the truncated Gross-Pitaevskii equation 
reached (microcanonical) equilibrium and that a conden- 
sation transition of the equilibrium was obtained when 
the initial energy was varied. 

The same condensation transition was later studied 
by Connaughton ct al. [2 7) and interpreted as a con- 
densation of classical nonlinear waves. Using a modified 
wave turbulence theory with ultraviolet cutoff, they ar- 
gued that the transition to condensation should be sub- 
critical. They found their theory in quantitative agree- 
ment with numerical integration of the GPE, using the 
same stochastic initial conditions than those of reference 
[13] . However, the authors later argued that, as weak 
turbulence theory is expected to breakdown nearby the 
transition to condensation, the subcritical nature of the 
transition predicted by their theory was not physical I28J . 

Berloff and Svistunov [53], starting from periodic ini- 
tial conditions similar to those of Davis et al. [15] . 
used a finite-difference scheme (exactly conserving en- 
ergy and particle number) to characterized the dynamical 
scenario of the relaxation toward equilibrium. Using the 



same finite-difference scheme, Berloff and Youd [5U] then 
studied the dissipative dynamics of superfluid vortices at 
nonzero temperatures and observed a contraction of the 
vortex rings that followed a universal decay law. 

Our main results are the followings. The classical abso- 
lute equilibrium of ideal fluids when generalized to GPE 
superfluids describes a standard [5TJ [55] second-order 
phase transition. Long transient with energy cascade 
and partial small-scales thermalization are present in the 
relaxation dynamics. Dynamical counter-flow effects on 
vortex evolution are naturally present in the system and 
the vortex ring have anomalous velocities caused by ther- 
mally excited Kelvin waves. 

The paper is organized as follows: Section[ll]is devoted 
to the basic theoretical background that is needed to ac- 
count for the dynamics and thermalization of the Fourier 
truncated GPE. 

In Sec. |III[ the thermodynamic equilibrium is ex- 
plored. The microcanonical and grand canonical distri- 
butions are numerically shown to be equivalent. Exact 
analytical expressions for the low-temperature thermody- 
namic functions are obtained. A standard second-order A 
phase transition is exhibited at finite-temperature using 
the SGLE-generated grand canonical states. 

In Sec. |IV[ the direct energy cascade is considered as 
a new mechanism for GPE thermalization. Using initial 
data with mass and energy distributed at large scales, a 
long transient with partial thermalization at small-scales 
is characterized. Vortex annihilation is observed to take 
place and is related to mutual friction effects. A bottle- 
neck producing spontaneous self truncation with partial 
thermalization and a time-evolving effective truncation 
wavenumber is characterized in the limit of large disper- 
sive effects at the maximum wavenumber of the simula- 
tion. 

In Sec. [V] the new SGLE algorithm is used to generate 
counter-flow states, with non-zero values of momentum, 
that are shown to be metastable under SGLE evolution. 
The spontaneous nucleation of vortex ring and the cor- 
responding Arrhenius law are characterized. Dynamical 
counter-flow effects are investigated using vortex rings 
and straight vortex lines arranged in crystal-like patterns. 
An anomalous translational velocity of vortex ring is ex- 
hibited and is quantitatively related to the effect of ther- 
mally excited finite-amplitude Kelvin waves. Orders of 
magnitude are estimated for the corresponding effects in 
weakly interacting Bose-Einstein condensates and super- 
fluid 4 He. 

Section IVII is our conclusion. The numerical methods 
and low-temperature thermodynamic functions are de- 
scribed in an appendix. 



II. THEORETICAL BACKGROUND 

This section deals with basic facts needed to under- 
stand the dynamics and thermalization of the Fourier 
truncated GPE. We first recall in section III A II the 
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(untruncated) GPE dynamics, its associated conserved 
quantities and the corresponding spectra; this material 
can be skipped by the reader already familiar with the 
GPE model of superflow [H [3] . The Fourier truncated 
GPE, its thcrmodynamical limit and the different statis- 
tical ensembles are then defined. 

The thermodynamics of the truncated system is intro- 
duced in section |II B| using the microcanonical distribu- 
tion. The canonical and grand canonical distributions are 
also used as they allow to directly label the equilibrium 
states by temperature and particle numbers. 

A stochastically forced Ginzburg-Landau equation 
(SGLE) is considered in section |II C and shown to de- 
fine a new algorithm that directly generates the grand 
canonical distributions. 



A. Galerkin truncated Gross-Pitaevskii equation 

1. Conservation laws and Galilean invariance of the GPE 

Superfhrids and Bose-Einstein condensates [3, 33 can 
be described at low temperature by the Gross-Pitaevskii 
equation (GPE) that is a partial differential equation 
(PDE) for the complex field ip that reads 



in 



where \ip\ 2 is the number of particles per unit volume, 
m is the mass of the condensed particles and g = , 
with a the s-wave scattering length. This equation con- 
serves the Hamiltonian H, the total number of particles 
N and the momentum P defined in volume V by 



H = 
N = 
P = 



d 3 a 



2m 



.9 



\ipY<Tx 



(tpVtp - ifiVi/j) d 3 x. 



(2) 
(3) 
(4) 



It will be useful for the next sections to explicitly write 
the conservation law of the momentum dt^(ipdjip — 
ipdjip) + d}Jlkj — 0, where the momentum flux tensor 
Ukj is defined, following ref. [6], as 



n 



ti 2 - - g 

i: i = ;r— ( d kipdjip + d k ipdjip) + S kj ( - 1 



4m 



v 2 H 2 



(5) 

It is well known that the GPE can be mapped into 
hydrodynamics equations of motion for a compressible ir- 
rotational fluids using the Madelung transformation de- 
fined by 



where p(x, t) is the fluid density and </>(x, t) is the ve- 
locity potential such that v = V0. The Madelung 
transformation ([6]) is singular on the zeros of if). As 
two conditions are required (both real and imaginary 
part of ip must vanish) these singularities generally take 
place on points in two-dimension and on curves in three- 
dimensions. The Onsager-Feynman quantum of velocity 
circulation around vortex lines ip = is given by h/m. 

When Eq.Q is linearized around a constant ^ = Aq, 
the sound velocity is given by c = g\A$\ 2 /m with dis- 
persive effects taking place for length scales smaller than 
the coherence length defined by 

£ = ^h 2 /2m\A \ 2 g. (7) 

£ is also the length scale of the vortex core [3j |6] . 

Following reference [5] we define the total energy per 
unit volume e to t = (H — fiN)/ V — [i 2 /2g where /i is the 
chemical potential (see section II B). Using the hydro- 
dynamical variables, etot can we written as the sum of 
three parts: the kinetic energy efci„, the internal energy 
ei n t and the quantum energy e q defined by 



eint 



1 

V 
1 

V 
1 

V 
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(8) 
(9) 
(10) 



Using Parseval's theorem, one can define corresponding 
energy spectra: e.g. the kinetic energy spectrum ekin(fc) 
is defined as the sum over the angles 



eidn(fc) 



<f i 



c v^ v 



k 2 dn h 



(ii) 



where dQ, k is the solid angle element on the sphere. The 
energy ekin can be further decomposed into a compress- 
ible part e£ in and an incompressible part ej^ by mak- 
ing use of the relation ^/pv = (^/pv) + (^/pv) 1 with 

V • (VP V )' = ( see M for details). 

Finally note that the GPE |T]) is invariant under the 
Galilean transformation 



?//(x, t) = i/;(x — vq^, t) exp 



im 
IT 



VG 



1 2 ' 

x v n t 

2 G 



Under this transformation Eqs.(2]|4) transform as 

H' = ]^mNv 2 G + P • v G + H 
N' = N 

P' = mNvG + P- 
2. Definition of the Fourier truncated GPE 



(12) 

(13) 

(14) 
(15) 



exp [i — < 

m n 



(*>*)] 



(6) 



For a periodical 3D system of volume V the Fourier 
truncated GPE is defined by performing a Galerkin trun- 
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cation that consists in keeping only the Fourier modes 
with wavenumbers smaller than a UV cut-off fc max . 
Expressing ip in terms of the Fourier modes as 



V>(x,t) = ^k(i)e lk x , with 



h ■ 

''■ 11111 



G If 



(16) 



and where /c m i n = 2~n jV 1 ^ is the smallest wavenumber. 
The Galerkin (Fourier) truncated Gross-Pitaevskii equa- 
tion (TGPE) is defined as 



..dA k h 2 k 2 v - , 

ki,k 2 



at 



2m 



where the Fourier modes satisfy A^ = if k > /c max 
and the sum is performed over all wavenumbers satis- 
fying |ki|,|k 2 |,|k 2 + ki|,|k + k 2 | < k max . This time- 
reversible finite system of ordinary differential equa- 
tions with a large number of degree of freedom M ~ 
(femax/fcmm) also conserves the energy, number of parti- 
cles and momentum. 

The direct numerical evolution of the convolution in 
Eq.( |17[ ) would be very expensive in computational time 
0{N^), where N is the resolution. This difficulty is 
avoided by using pseudo-spectral methods |34| and the 
non-linear term is calculated in physical space, using 
FFTs that reduce the CPU time to <3(A 3 logA). Intro- 
ducing the Galerkin projector Vo that reads in Fourier 
space pG^k] = 0(k max — k)A-^ with #(•) the Heavside 
function, the TGPE |l7|) can be written as 



ih 



dip 

at 



■V G [ 



2m 



V 2 V + 57>g[M 2 ]# (18) 



Equation (181 exactly conserves energy and mass and, if 



it is correctly de-aliased using the 2/3-rule [33] (dealias- 



ing at k n 



2N 
3 2 



), it also conserves momentum 



(see Appendix [A] for a explicit demonstration). The 
Galerkin truncation also preserves the Hamiltonian struc- 
ture with the truncated Hamiltonian given by H = 



Let us remark that perhaps a more standard defini- 
tion of dealiasing in Eq.|l8]) could have been "Pg[|V , | 2 '0] 



using 1/2-rule (dealiasing at fe max = |^r) rather than 
V G [P G [\ip\ 2 }Tp] with the 2/3-rule. Using the former defi- 
nition removes the restriction |k 2 | < fc max on the con- 
volution in Eq.(17|. Both methods are equivalent in 



the partial differential equation (PDE) limit (exponen- 
tial decay of energy spectrum for k <C fc max ) and admit 
the same invariants. However the scheme of Eq.(18) is 
preferable because fc max is larger at the same resolution. 
If dealiasing is not preformed in equation ( 18 ) the errors 



in the conservation of momentum can rise up to 50% in 
a few units of time (see Appendix [A|. In a finite differ- 
ence scheme the conservation of momentum should also 
be checked carefully as it is bound to produce spurious 
effects. 



Another effect caused by periodic boundary condition 
is that the velocity vq in the Galilean transformation 
([121 is quantized by the relation 



H 2ir 

"7/1/3 



no, 



(19) 



where nc G Z 3 and v s becomes continuous only in the 
limit h/ (mV 1 / 31 ) — > 0. The Galilean invariance is slightly 
broken by the TGPE (17) because of modes close to the 
truncation wavenumber fc max . However it is recovered 
in the PDE limit where high wavenumber modes are 
converging exponentially and also in the thermodynamic 
limit: — » oo defined below because the offending 
terms represent only a surface effect in Fourier space. 



3. Thermodynamical limit and statistical ensembles 

Let us first notice that the energy H, the number of 
particles N and the momentum P in Eqs.(2 4) are all 
proportional to the total number of modes A 



and therefore are all extensive quantities. Also note that 
by definition of the coherence length Q, the number 
£fc max determines the amount of dispersion at truncation 
wavenumber in the system. 

The thermodynamic limit V — > oo of the truncated 
Gross-Pitaevskii system is thus defined as the limit 



£fcm ax = constant, 



(20) 



in order to obtain equivalent systems. In this limit the 
relevant thermodynamic variables are the intensive quan- 
tities H/V, N/V and P/V. In practice, to perform nu- 
merical computations we will fix the volume to V = (27r) 3 
and we will vary fc max (see paragraphs before section III I . 

Let us define, as usual the microcanonical ensemble 
[35) by the probability dw of finding the system in states 
with given values of energy H m , number of particles N- m 
(the subscript "in" stands for initial data) and momen- 
tum P; n given by: 

dw = constant e s 5{H~H in )5(N-N in )6 3 (P-P in )dHdNd 3 P, 

( 21 ) 

where S = log T is the entropy with T the number of 
accessible micro-states. 

Microcanonical statistical states can be obtained nu- 
merically by time-integrating the TGPE until the sys- 
tem reaches thermodynamic equilibrium |T3l El] . These 
thermalized states are formally determined by the con- 
trol values .ffin, A; n and Pj n that are set in the initial 
condition. It has been shown in references [13 [57] by 
varying the values of H- m that TGPE present a phase 
transition analogous to the one of Bose-Einstein conden- 
sation, where the amplitude at 0-wave-number Aq vanish 
for finite values of H lri . Let us remark that an explicit 
expression of dw or S cannot be easily obtained in the 
microcanonical ensemble and therefore the temperature 
is not easily accessible. 
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A simple way to explicitly control the temperature 
is to use the canonical or grand canonical formulation. 
The grand canonical distribution probability is given by 
a Boltzman weight 



F 



1 

— < 
Z 

H 



-/3F 



fxN-W-P, 



(22) 
(23) 



where Z is the grand partition function, j3 is the inverse 
temperature and [i is the chemical potential. In what 
follows we will refer to W as the counterflow velocity. 
Note that when W = 0, F = H — fiN and the statistic 



weight of distribution (22) corresponds to the so-called 
\(j> 4 theory studied in second order phase transitions [3U 
132] . This point will be further discussed in subsection 

mfci 

Finally remark that the states with W ^ are ob- 
tained, in the thermodynamic limit, by a Galilean trans- 
formation of the basic W = state (see below Eq.(72|). 
However, for finite size systems, because of the quantifi- 
cation of the Galilean transformation (Eqs.(12) and (19)) 
new metastable states with counterflow appear. These 
metastable states and their interactions with vortices will 
be studied in detail below in section IV Al 

In the grand canonical ensemble ( 22p3 ) the mean en- 
ergy H, number of particles N and momentum P are 
easily obtained by defining the grand canonical potential 



and using the relations 



N = - 



on 



on 
&p : 



H = 



w 



(24) 



-HN + W-P. (25) 



Observe that the microcanonical states (21) are char- 
acterized by the values H m , N- ln and Pi n . On the other 
hand, the grand canonical states are controlled by the 
conjugate variables: j3, fi and W. The different sta- 
tistical ensembles are expected to be equivalent in the 



thermodynamics limit ( 20 ) and therefore 



Hhi = H , N ia = N, P in =P, 



(26) 



in this limit. The equivalence of ensembles will be nu- 
merically tested below in subsection |III A| 

In the grand canonical ensemble, the pressure p is usu- 
ally defined from the grand canonical potential ( 24 ) by 
the relation [35] f2 = — pV. This definition presents two 
problems in the TGPE system. First, due to classical 
statistics has a logarithmic divergence at (3 = oo. Sec- 
ond, this definition does not coincide with the standard 
relation in fluid dynamics involving the diagonal part of 
the momentum flux tensor Ily (see Eq.(5|). Both these 
problems can be solved by considering the total number 
of modes as a new thermodynamics variable, as we will 
sec in the next section. 



B. Thermodynamics of the truncated system 

When a Galerkin truncation is performed on a system 
a new variable k ma _ x explicitly appears. One thus find 
that the thermodynamic potentials depend on the total 
number of modes. Denoting Xj\f the conjugate variable 
to the total number of modes Af the standard thermody- 
namic relation for the energy easily generalizes as 

dE = -pdV + TdS + pdN + X^dM + W -dP (27) 

with S the entropy and where we have included the to- 
tal momentum dependence dP. As in Landau two-fluid 



model lj Eq. ( 27 1 is written in a system of reference where 
v s = V(f) = 0~[the bar standing for some ensemble av- 
erage) and E = H is the macroscopic energy. |67j . We 
will omit the bar over the others microscopic quantities. 
Note that the Fourier modes formally play the role of 
"particles" and X_\f is formally the "chemical potential" 
associated to those "particles" . 

The thermodynamic potentials can be easily general- 
ized to take in to account the new variables. It is useful 
to define the Gibbs potential G, grand canonical il and 
a generalized grand canonical potential f2' (with a Leg- 
endre transformation on Af) as 

G = E - TS + P V - W-P (28) 
n = E - TS - fiN - W-P (29) 
fl' = E - TS - fiN - Aa/A/" — W-P (30) 

from where their respective variations follows: 

dG = Vdp - SdT + p,dN + X^/dAf — P ■ dW (31) 
dfl = -pdV - SdT - Ndfi + Xj^dN - P ■ dW(32) 
dQ! = -pdV - SdT - Ndfi - NdX M - P ■ dW{33) 

Based on standard arguments of extensive variables 
[3 5) and noting that Aw and W are intensive variables 
we find the standard formula of the Gibbs potential with 
two types of particles 



G = fiN + X^JV. 



(34) 



Using Eqs.((28|) and |34j) in Eqs.(|29| and (|30j) we find 

n = - P v + x^m, n' = -pv (35) 



The relations ( 27||35 1 determine all the thermodynamic 
variables and potentials. For instance the pressure p can 



be obtained from Eq.(32), Eq.(33l or Eq.(35) by 



dV 



7>,AA,W 



tt-X M N 
V 



9l 

V 



(36) 



where X M = j§\ v ^ w - 

We proceed now to show that thermodynamic defini- 
tion ( 36 1 of the pressure coincides with the standard re- 



lation in fluid dynamics. In order to make explicit the 
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dependence of the energy H on the volume V let us de- 
fine the dimensionless space variables x = x/Vz and tp = 
V x / 2 ifj. Expressed in term of these variables the Hamil- 

tonian (2) reads H = J d 3 x (jjt A-|V^| 2 + if |^| 4 ) 

Taking the derivative with respect to V and reintroduc- 
ing x and ip yields 



dH 

W 



v 1 d " x 



h 2 2 

2m 3 



(37) 



This expression corresponds to the spatial average of the 
the diagonal part of ILfc (see Eq.Q). As by definition 
E = H and the derivative has been implicitly done at 
constant total number of modes and momentum we find, 



using the thermodynamic relation (27 1 and Eq.(37), that 
the pressure satisfies 



dE 
dV 



S,N,Af,P 



dH 
dV 



iVJV.P 



(38) 



where the second equality holds for adiabatic compres- 
sions 1551. 



Finally by replacing ft in Eq. ( 29 1 we obtain the ther- 
modynamic relation 



E + pV - fxN - W P = TS + XjvAf. 



(39) 



Let us remark that, in a classical system, the entropy is 
defined up to an additive constant related to the nor- 
malization of the phase-space. However the quantity 
TS+Xj^Af is completely determined because each term in 
the left hand side of Eq. ( 39 ) is well defined. By the same 



arguments d (AfXj\f /T) is also a completely determined 
quantity. If the variable Af had not been taken into ac- 
count, the corresponding pressure would be —Cl/V and 
therefore wrongly defined and depending on the normal- 
ization constant. The grand canonical potential f2 will 
be explicitly obtained at low-temperature in subsection 
|IIIB| where the above considerations can be explicitly 
checked. 



C. Generation of grand canonical distribution 
using a stochastic Ginzburg-Landau equation 

Grand canonical equilibrium states are given by the 
statistics ( 22p3 ). They cannot be easily obtained be- 
cause the Hamiltonian H in Eq.([2| is not quadratic 
and therefore the statistical distribution is not Gaussian. 
Nevertheless it is possible to construct a stochastic pro- 
cess that converges to a stationary solution with equilib- 
rium distribution ( 22p3| . This process is defined by a 
Langevin equation consisting of a stochastic Ginbzurg- 
Landau equation (SGLE) that reads 



~W 

(C(x,t)CV,f)) 




<5(i-*')<5(x-x'), 



C(k,t) (40) 
(41) 



where F is defined in Eq.(23l and ((k, t) is the [k 



truncated) Fourier transform of the gaussian white-noise 
C(x, t) defined by Eq.(41). The Langevin equation ( [40] - 
pdj) explicitly reads in physical space 



<9V 



h 2 

^g[^V 2 V + ^ - gV G [ 

2m 



2 ]V> - ihw ■ W] 




Pg[C(x,*)]. 



(42) 



In the P — > oo limit Eq.(42) reduces to the advective 



real Ginzurg-Landau equation (up to a redefinition of /i) 
that was introduced in reference [5]. This equation has 
the same stationary solutions of than the TGPE ( 18 ) in 



a system of reference moving with velocity W. When the 
term is also included in the TGPE it has, because of 
particle number conservation, the only effect of adding a 
global time-dependent phase factor to the solution. 

The probability distribution P [{A k ,_A? }k<k max ] of the 
stochastic process defined by Eqs.( 40]|41 1 can be shown 
to obey the following Fokker-Planck equation [35J |37] 



at 



E 



k<k n 



d 

8A* 



1 OF 
VhdAt 



1 <9P 
VhpdA^ 



+ c.c. (43) 



It is straightforward to demonstrate that the probabil- 
ity distribution ( 22 ) is a stationary solution of Eq. ( 43 1 , 
provided that pF is a positive defined function of 
{^k,^4 k }k<k max - 

If one wishes to directly control, instead of the chemical 
potential /z, the value of the number of particles N or 
the pressure p, the SGLE must we supplied with one of 
two ad-hoc equation for the chemical potential. These 
equation simply read 



d/i 

dt 
d/i 

dt 



-u N (N-N*)/V 
-v p {p-p*) 



(44) 
(45) 



where the pressure p is computed as p ~ (see 
Eq.(37)). Equation (44) controls the number of particles 
and fixes its mean value to the control value N* . Simi- 



larly Eq.(|45|) controls the pressure and fixes its value at 
45 1 are not compatible and they must 
Depending on the type of 



p* . Equations (44 
not be used simu 



taneously. 

temperature scans, the SGLE must be used together with 
either Eq.(44), Eq.(45) or alone with a fixed value of /.i. 



In the rest of this paper we will perform several numeri- 
cal simulations of the TGPE (IT7|) and SGLE |42l. For nu- 



merics, the parameters in SGLE (omitting the Galerkin 
projector Vq) will be rewritten as 



dt 



= avV 2 ip + VLQip- P^ii-iW ■ Vip- 



k B T 
a 



-c, 



with similar changes for TGPE. 

In terms of ao, Slo an d Po the physical relevant pa- 
rameters are the coherence length £ and the velocity of 
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sound c denned in section II A 1 (see Eq. ^ and the text 
before). They can be are expressed as 



V 2a oPof 



(46) 



with p* = Oo//?o- The value of the density at T = 
set to p* = 1 in all the simulations presented below. In 
order to keep the value of intensive variables constant 
in the thermodynamic limit ( 20 1 , with V constant and 
J; max — ► oo the inverse temperature is expressed as = 
l/ktfT where fc/y = V/N. With these definitions the 
temperature T has units of energy per volume and 47rao 
is the quantum of circulation. 

With £ fixed, the value of £/c only determine a time- 
scale. The velocity of sound is (arbitrarily) set to c = 2 
and the different runs presented below are obtained by 
varying only the coherence length £, the temperature T, 
the counterflow velocity W and the UV cut-off wavenum- 
ber fc ma x- The number £fc max is kept constant (at the 
value £fc m ax = 1-48) when the resolution is changed, ex- 
cept in section [TV B| where dispersive effects are studied 
(using a larger £fc max ). This choice of £ ensures that vor- 
tices are well resolved (e.g. compare Figjfjja below with 
Fig. 12 of ref. 6J). In the present work we use resolutions 
varying from 32 3 to 512 3 colocation points (fc max = 10 
to 170 respectively). Finally in all numerical results the 
energy and momentum are presented per unit of volume 
V = (27r) 3 and the control values of number of particles 
and pressure in Eqs.( 44p5 1 are set to mN* /V = p* = 1 
and p* = c 2 p* 2 /2 = 2. Numerical integrations are per- 
formed with periodic pseudo-spectral codes and the time- 
stepping schemes are Runge-Kutta of order 4 for TGPE 
and implicit Euler for SGLE. 



III. CHARACTERIZATION OF 
THERMODYNAMIC EQUILIBRIUM 

In this section, the thermodynamic equilibrium is ex- 
plored and characterized. The microcanonical and grand 
canonical distributions are first shown to be numerically 
equivalent in a range of temperatures by comparing the 
statistics of GPE and SGLE generated states in section 
|III A| The steepest descent method is then applied to the 
grand partition function in section |III B| to obtain exact 
analytical expressions for the low-temperature thermody- 
namic functions. The basic numerical tools are validated 
by reproducing these low-temperature results. In section 
|III C| a standard finite-temperature second-order A phase 
transition is exhibited using the SGLE-generated grand 
canonical states and the deviations to low-temperature 
equipartition are characterized. 



Comparison of microcanonical and grand 
canonical states 



We now numerically compare the statistics of the grand 
canonical states produced by the new algorithm SGLE 



to the statistics of the microcanonical states obtained by 
long-time integrations of TGPE. The coherence length 
is set to £ = y/2/10 and 32 3 collocation points are used 
(&max = 10). The initial condition for the TGPE runs 
are chosen with random phases in a similar way than 
in references [l3|[27]. We obtain low, medium and high 
values of the energy with constant density p — mN /V — 
1 (see table [I). 

TABLE I: Parameters of TGPE initial condition and time 
steps. 



H 


T 


TGPE time steps 


SGLE time steps 


0.09 


0.09 


40000 


9600 


0.5 


0.5 


20000 


9600 


1.96 


1.8 


20000 


9600 


4.68 


4 


20000 


5000 



To compare with the SGLE generated statistics a scan 
in temperature at constant density p = 1 is performed 
in order to obtain the temperature corresponding to the 
energies of the TGPE runs. Using the thermalized final 
states obtained from TGPE and converged final states of 
SGLE histograms of the of the density p(x) in physical 
space are confronted in Figjl] They are found to be in 
excellent agreement. 
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FIG. 1: (Color online) Comparison of density histograms ob- 
tained by SGLE and TGPE dynamics (£ = 2/Wy/2 and res- 
olution 32 3 ) with energy equal to a) H = 0.09, b) H = 0.51, 
c) H = 1.96 and d) H = 4.68 (see table [j}. The solid line in 
a) is a Gaussian of standard deviation 5p? = 0.016 (see below 
Eq.(60l) computed with the low-temperature calculations of 
section IIIIBI 



Observe that when the energy (or temperature) in- 
creases more weight becomes apparent on the histograms 
near p = 0, indicating the presence of vortices. The 



Gaussian character of the histogram in FigjT] a motivates 
the low-temperature calculation of the next section. Ob- 
serve that, even at this relatively low 32 3 resolution, the 
thermodynamic limit has been reached in the sense that 
the micro and grand canonical distribution coincide. We 
can thus safely use, at 32 3 and higher resolutions, the 
SGLE to prepare absolute equilibria of the TGPE. 

Let us finally remark that the SGLE numerically con- 
verges much more rapidly toward absolute equilibrium 
than the TGPE, as displayed on table [TJ Also taking into 
account the (computationally expansive) accurate con- 
servative temporal scheme needed for the integration of 
the TGPE, the SGLE yields a (very) large economy of 
the CPU time needed to reach equilibrium. On the lo- 
cal machines where these computations were performed 
the SGLE was typically more than 10 times faster than 
TGPE. 



p = Kk, at leading order F can be rewritten as F 
F +Fi + F 2 with 



r9 ]A |4 
P 2 



Fo = F(||Ao| 4 -/i|A | 2 ) 



(49) 



F i = ^]T(^-/i + 2 ff L4o| 2 -W.p)|A p | 2 (50) 
* — ' 2m 



p#0 

2 



F 2 = V^J2 A * 2 o A p A -p + A 1 A 1 A * 

p#0 



(51) 



In order to obtain the low-temperature partition func- 
tion we need to compute the determinant of the matrix 
d ® q A — /xo^(5 Pl o<5q,o- This determinant can be obtained 
by making use of the Bogoliubov transformation 



.4. 



u p B p 



(52) 



B. Low-temperature calculation 

The gaussian histogram of FigjTJa strongly suggest 
that some quadratic approximation should be able to ob- 
tain exact analytical expressions for the thermodynamic 
functions at low temperature. In this section we use 
a such an approximation to compute the grand parti- 
tion function Z and the grand canonical potential [68] 
= -/3" 1 log.Z defined in p9b. 



The first step is to express the energy F of Eq.(23) 
in terms of the Fourier amplitudes A^. This leads to 
a non quadratic function F [At, A^] explicitly given in 
appendix [B] (Eqs |Bl|B3[ ). The grand partition function 
is a product integral over all the Fourier amplitudes 



with u p = r^^jf, v p = Ja^^I and where L p 

is determined by imposing the diagonalization of F — 
Ho\Ao\ 2 V. L p is explicitly given in Eq.(B6|. It is easy 
to show that (52 1 is a canonical transformation and that 



the normalization condition of the corresponding Poisson 
bracket implies \u p \ 2 — \v p \ 2 = 1. 

Expressing F in the Bogoliubov basis we obtain 



F = V 



||A | 4 - Ml^ol 2 + £ (c(p; A*,A«o) - W • p) |i? p | ; 

p#0 

(53) 

with the dispersion relation (see appendix |B[) 



dA dA* 
2^ 



TT dA ^ dA *^ c -0F\A k ,At] e(p;/i,Mo) - \j ( M + 2/i + 



k<k„ 



2tt 



2m 



(/i + Mo) 2 . (54) 



The integrals in (471 cannot be done explicitly, however 



it is possible to give a low-temperature approximation us- 
ing the method of steepest descent [31J [38] . We also add 
to F an external field with value — /io|Ao| 2 ^ in order to 
explicitly obtain the mean value of condensate Fourier 
mode \Ao\ 2 by direct differentiation. The physical parti- 
tion function is finally obtained by setting [1q = 0. The 
integrals are dominated by the saddle-point determined 
by J^pr — /iqAqVS^o = that yields the solution (see 



Let us recall that the exited modes B p are called phonons 
in quantum mechanics. In the present classical case, be- 
cause of classical statistics and quadratic Hamiltonian, 
there will be equipartition among phonon modes. Re- 
placing the value of the chemical potential by its saddle- 



point expression fi = g\A \-\ (at no = 0), Eq.(54) 



sp 



Eqs.(TB4| and (TB5J) 



9\Ao\ 



fJ> + Mo 



-4, 



for k ^ 0, (48) 



where the subscript "sp" stands for saddle-point. Note 
that in general \Aq\ 2 ^ |j4o| 2 | and the mean value is 

Isp 

equal to the saddle-point one only at T = 0. Other solu- 
tions that can be obtained when W ^ will be discussed 
in detail in section [V] 

At the saddle-point (48 1 A^ vanishes for k ^ 0. We 



thus need to keep only quadratic terms in A^ to obtain 
the low-temperature approximation. Using the notation 



yields the (standard, see ref.[39]) Bogoliubov dispersion 

relation e(p) = P\J 9 ^^ - + 4—2- Note that e(p) can also 
be directly obtained from the GPE by expressing ip in hy- 
drodynamics variables, using the Madelung transforma- 
tion ^ and linearizing around an homogenous density 

The partition function now trivially factor- 
izes in independent parts Z(/3, V, /1, W, A/", /io) = 
Z (/3, /j,, /i ) Ilp^o Zp(P> M> w , Mo) where 



Z (/3,V,n,Ho) = V2^J—e ^ 




K e iPl M,Mo) - W • p) 



(56) 
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The total number of modes Af = Ylk 1 anc ^ the g ran d 
canonical potential 



thermodynamic relation ( 32 I 



logZ + ^2 logZp 

p#0 



Replacing the sum by an integral the expression for 
the number of modes reads 



(57) 



TV 



P 2 V 
2ir 2 h 3 



dp 



P 3 V 

max v 



6n 2 h 3 



(58) 



are sums over all wave-numbers from which all thermo- 
dynamic quantities can be directly obtained by using the 



Setting W = (0,0, w) the integral form of Eq.(57) reads 



V{H + Mo) 2 

V(p + 2p Q ) 
2 5 



(m + Mo) 2 - Ppwz 



-1 

pL x v { 2 

6tt 2 /3/i 3 \3 



log [j3e(P max ;t,)}-f 



Amp 

max 



2/t 



Mo 



fo 



dz dp 
2 

Amp 

max _ 



(59) 



In order to obtain E q.(|59| ) the thermodynamic limit ( 20 ) 
of infinite volume [69 was taken and the conditions 
p/m, Mo/m ^ 1 were used. The functions f[z] 



and fo[z] are explicited in Eqs.( Bl6"]|Bll ). Note that 
the dependence of the grand canonical potential Q on 
the number of modes Af is implicitly given by P max and 
Eq.(58|. The first term in 17 is due to the condensed 
mode at p = 0. 

The low-temperature approximation to all thermody- 



namic functions is directly obtained from equation (59) 



by first setting p — and then differentiating ( 59 ) , using 



relation ( 32 ) . It is straightforward to check that both def- 



inition of the pressure in Eq.(36) coincide. Furthermore 



the higher order moments of the density can be easily 
computed by taking successive derivatives of the grand 
canonical potential. For instance it is straightforward to 
show that the variance of the density p (see solid line on 
histograms displayed on Fig{l]a) is given by 



V 2 (5 P 2 ) = -13- 



, d 2 n 

dp 2 



(60) 



It can also be checked on the explicit expression for the 



entropy S (see Eqs.(B9l) that, as expected for a classical 



system, the entropy depends by a logarithmic term on the 
phase-space normalization. However the function TS + 
\_\fj\f is independent of phase-space normalization (see 
discussion below Eq.([39|). 

Finally, low-temperature expressions for the energies 
(8p0 ) and their corresponding spectra can be easily ob- 
tained using Madelung's transformation ^ . At low tem- 
peratures the fluctuations are smalls and e^in depends 
only on </> and e q + e- m t only on p. The total energy is 
thus decomposed in two non-interacting terms. Equipar- 
tition of energy between the total kinetic energy ekm and 
quantum plus internal energy e q + e mt is thus expected 
at low temperature. 



The next subsections will be concerned with the van- 
ishing counterflow case w — 0. The states with non-zero 
counterflow w will be studied in details in section Ivl 



C. A transition and vortices 

To characterize the condensation transition, we present 
here four temperature scans performed using SGLE ( |42| . 
Three of them are at resolution of 64 3 with respectively 
constant chemical potential, density and pressure (us- 
ing Eqs.( 44p5] )). The fourth scan is performed at con- 
stant pressure but at a resolution of 128 3 . The coherence 
length is fixed so that £fc max = 1.48 is kept constant. 

Figure [2]a displays the results of the scans. Observe 
that the low-temperature behavior is in good agreement 



with the analytical calculations of section IIIB and the 
explicit formulae given in appendix [B] Also observe that 
the constant pressure scans at resolutions of 64 3 and 128 2 
coincide for all temperatures showing that the thermody- 
namic limit ( [20] ) discussed in section II B is obtained at 
these resolutions. 

Figure [2jb displays the temperatures dependence of 
the condensate fraction |^4o|/p for the four SGLE runs. 
Observe that the condensation transition previously ob- 
tained (in the constant density case) by microcanonical 
simulations in references [T31 [23 i s reproduced and also 
present (at different critical temperatures) in the con- 
stant pressure and chemical potential scans. 

The SGLE algorithm directly provides the tempera- 
ture as a control variable. It thus allows to easily obtain 
the specific heat from the data. Figure [2|c displays the 
specific heat at constant pressure c p — ^ | for the scan 

at resolution 128 3 . Let us remark that the (w = 0) statis- 
tic weight of distribution ( 22|23 1 corresponds to that 
of (standard two-component) second order phase tran- 
sitions [3TJ [31] • We thus expect the condensation transi- 
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FIG. 2: (Color online) a) Temperature dependence of the den- 
sity p, pressure p and chemical potential p for SGLE scans 
at constant density, pressure and chemical potential (see leg- 
end on figure), b) Temperature dependence of the conden- 
sate fraction \Ao\ 2 /p (same scans as in a)), c) Specific heat 
at constant pressure and resolution 128 3 the solid 



on 

8T 



line corresponds to a fit (see Eq.(61 1). d) Temperature depen- 
dence of the energies e£ in , e^, ekin and e q + ei nt at constant 



density; equipartition of energy between ekh 
apparent at low temperatures. 



and e q + eint is 



tion visible on Fig[2]c to be in this standard class. This 
point is confirmed by the solid lines in Figure [2]c that 
correspond to a fit with the theoretical prediction given 
by the renormalization group (RG) 



«RG 



: (1 



-a±M A 



_|2A 



B ± (61) 



where r = T j? x and the + and — signs refer to T > T\ 
and T < T\, see reference [30]. The fit was obtained in 
the following way: first the identification of the transi- 
tion temperature T\ was done by finding the zero of the 
linear interpolation of the second order difference of H, 
discarding the three closest point to the zero of |A | 2 //°- 
Then, using the critical exponents cvrg = —0.01126 and 
A = 0.529 given by the RG, the data was fitted as in 
reference |40j over the non-universal constant. The ob- 
tained values are A + /A~ — 1.42, to be compared with 
1.05 in reference |40j . This discrepancy is probably due 
to finite size effects. 

Finally on Fig[2]d the temperature dependence at con- 
stant density of the different energies (JsJJTO ) expressed in 
terms of hydrodynamical variables is displayed. Observe 
that the incompressible kinetic energy EL vanishes for 
low temperatures T <C TC , where Tf = 3.31 is the tran- 



sition temperature at constant density. This vanishing is 
connected to the disappearance of vortices, that it is also 
manifest in the density histograms in Fig[T] At low tem- 
perature equipartition of energy between the total kinetic 
energy e^n and quantum plus interna l ener gy e q + e lnt 
(see the discussion at the end of section III B ) is apparent 
on Figj2}d. 



IV. ENERGY CASCADE, PARTIAL 
THERMALIZATION AND VORTEX 
ANNIHILATION 

A new mechanism of thermalization through a direct 
cascade of energy is studied in section |IV A| Using ini- 
tial conditions with mass and energy distributed at large 
scales, a long transient with partial thermalization of the 
density waves is obtained at small-scales. Vortex annihi- 
lation is observed to take place and is related to mutual 
friction effects. A bottleneck effect that produces spon- 
taneous self truncation with partial thermalization and a 
time-evolving effective truncation wavenumber is charac- 



terized in section IV B for large dispersive effects at the 
maximum wavenumber of the simulation. 



A. Partial thermalization 

We now study the (partial) thermalization of the su- 
perfluid Taylor-Green (TG) vortex. This flow, that was 
first introduced in reference [BJ, develops from an initial 
condition that is prepared by a minimization procedure 
using the advected real Ginzburg-Landau equation (AR- 
GLE) [BJ. The nodal lines of the initial condition ipTG 
are the vortex lines of the standard TG vortex and obeys 
all its symmetries. Numerical integrations are performed 
with a symmetric pseudo-spectral code, making use of 
the TG symmetries to speed up the computations and 
optimize memory use, as described in reference [BJ. We 
use the equivalent to 64 3 , 128 3 , 256 3 and 512 3 collocation 
points and the coherence length is set such that, in all 
cases, £fc max = 1-48. 

Vortices and density fluctuations corresponding to the 
512 3 run are visualized on Fig{3] using the VAPOR [70] 
software. The short time behavior, see Fig[3]a-c, cor- 
responds to the GPE superfluid turbulent regime pre- 
viously studied in ref.[5J. A new TGPE thermalization 
regime where vortices first reconnect into simpler struc- 
tures and then decrease in size with the emergence of a 
thermal cloud is present at latter times, see Fig(3]d-e. 

To further study this new TGPE regime, the temporal 
evolution of e kin , e\. in , e£ in , e q + e^t is shown on Figjija 
and the corresponding energy spectra are displayed on 
FigJH] Observe that, at t = 0, ej^ contains almost all 
the energy because of the highly vortical initial condition. 

The early times (t < 15) behavior corresponds to the 
PDE regime of the GPE (nl) that was previously reported 
in references [SI [B] . An energy transfer is observed from 
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FIG. 3: (Color online) 3D visualization of density at t = 0, 5, 
10, 20, 31 and 55 at resolution 512 3 . Vortices are displayed as 
grey (red) isosurfaces and the grey (blue) clouds correspond 
to density fluctuations. 



e^. in to the other energies (e£ in and e ln t + e q ) that are 
associated to the density waves. 

Continuing the temporal integration the spectral con- 
vergence to the GP partial differential equation is lost. 
The dynamics becomes influenced by the truncation 
wavenumbers fc max and thermalization starts to take 
place. Two new regimes are observed. The first one for 
20 < t < 80 corresponds to a partial thermalization at 
small-scales: see Fig[5jb-d. Observe that equipartition of 
e£ in and e q -I- e int begins to be established in this phase. 
The thermalized zone then progressively extends to larger 



wavenumbers. During this phase eL decrease at almost 
constant rate (see Figj4|a) . As shown on Fig|4jb, this 
phase is delayed when the resolution is increased at con- 
stant £A; max . 

Around t — 80 (Figj4ja and[5jd) equipartition is estab- 
lished for each wave-number and e^ in almost vanishes. 
This vanishing is related to the disappearance of vortex 
lines that first reconnect into simpler structures which 
then decrease in size and number and finally disappear 




100 



100 



FIG. 4: (Color online) a) Temporal evolution of energies e^, 
e kim e kin an d e q + eint at resolution 256 3 . At large times, the 
incompressible energy vanishes and equipartition of energy 
between ekin and e q + eint is observed, b) Temporal evolution 
of e^ at resolution of 64 3 , 128 3 , 256 3 and 512 3 with constant 
£fcmax = 1.48. 



(as can be directly observed on the density visualizations 
corresponding to the 256 3 run, pictures not shown). Note 
that the annihilation of the vortices can also be related 
to the contraction of vortex rings due to mutual fric- 
tion reported in ref. 30J. For t > 80 the system finally 
reaches the thermodynamic equilibrium. The final ab- 
sence of vortices and equipartition of energy between e£ in 
and e q -I- e- m t , as can be directly checked on the temper- 
ature scan in Fig(2jd, is a consequence of the low energy 
of the initial condition "0tg- 

We have thus presented for the first time a new mecha- 
nism of thermalization through a direct cascade of energy 
of the TGPE similar to that of the incompressible trun- 
cated Euler equation reported in reference |21j . 



B. Dispersive slowdown of thermalization and 
bottleneck 



We now turn to the study of dispersion effects on the 
thermalization of the TGPE dynamics and on vortex an- 
nihilation. To wit, we prepare three different initial con- 
ditions with different values of £fc max using the TG ini- 
tial condition described in the preceding section. We 
fix the value of the coherence length to £ = y/2/20 and 
use resolutions of 64 3 , 128 3 and 256 3 corresponding to 
(fcmai = 1.48, 2.97 and 6.01 respectively. The three ini- 
tial condition therefore represent the same field at differ- 
ent resolutions. 

The temporal evolutions of ekin, ej (in , e£ ir and e Q 



-eint 
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for the three runs (indexed by the resolution) are dis- 
played on FigjBja. They are identical until t s» 5 where 
the run of resolution 64 3 starts to lose its spectral con- 
vergence. At at t rj 20 all runs appear to have ther- 
malized on Fig(6ja. However the kinetic energy spec- 
tra on Fig(6jb shows a clear difference between the runs 
(the dashed line corresponds to k 2 power-law scaling). 
The high-wavenumber modes of the 64 3 run are thermal- 
ized. For the 128 3 run the high-wavenumbers begin to 
fall down and, at resolution 256 3 , two zones are clearly 
distinguished. An intermediate thermalized range with 
an approximative k 2 power-law scaling is followed by a 
steep decay zone well before fc max — 85. Remark that in 
the 256 3 run the spectral convergence is still ensured and 
the (partial) thermalization is thus obtained within the 
GP PDE-dynamics. 

The temporal evolution of ekin(^) for the 256 3 run is 
displayed in Figj6]c. The large wave-number k~ 3 power- 
law behavior at t — is an artifact of the high-fc de- 
composition of energies in the presence of vortices (see 
pp. 2649-2650 of ref.[6 and [JT]) and a faster decay is 
recovered as soon as the vortices disappear. The ther- 
malized intermediate zone is observed to slowly extends 
to smaller wave-numbers. This naturally defines a self- 
truncation wave-number k c (t) where the energy spectrum 
starts to drastically decrease. 

In order to determine k c (t) we have tested fits to 



ekin(fc) using two type of trial spectra with three 
free parameters: eg t i(fe) = A(t)k~ n exp [— 26(t)k] and 
e fit Ti(fc) = A(t)k- n exp[-j(t)k 2 ]. The e fit n(fc) fit was 
found to work better in the sense that it both gives the 
correct n = — 2 prcfactor at intermediate and large times 
and also gives a better fit to the data at high k (data not 
shown). Fixing the prefactor at the value n = —2, we 
finally define our working two-parameter fit as: 



entile, t) = A{t)k 2 



' [( 16 ) 3 ( kjjt) ) ] 



efit(fcj t)dk. 



(62) 
(63) 



The factor (9tt/16) 1/3 in Eq.fl62| was set in order to 
obtain the limits Ak^/3 and ^4fc 3 lax /3 for eet(i) when 
fcmax ~ > oo and k c — > oo respectively. The fits are also 
displayed in Fig[6jc. They are in good agreement with the 
data after vortices have disappeared. The temporal evo- 
lution of eet(f) is displayed in Figj6]a. It does converge to 
the thermalized value of the energy. Finally the temporal 
evolution of the self-truncation wavenumber k c (t), which 
seems to have a well defined limit at infinite resolution, 
is displayed in Fig(6jd for the three runs. 

An open question is wether k c is bounded in time in 
the PDE regime where k c <C fc max - In other words, is 
thermalization of the £fe max 3* 1 truncated system sim- 
ply delayed or completely inhibited when £fc ma x is large 
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a) 









e kin 256 

ei. 256 3 



wavenumber is determined using the integral formula 




FIG. 6: (Color online) a) Temporal evolution of energies (as 
in Fi^5]a) for £fc max = 1.48, 2.97 and 6.01 (resolution 64 2 , 
128 3 and 256 3 respectively). Yellow stars are the kinetic en- 
ergy reconstructed from fit data using Eq.( |63[ |. b) Kinetic 
energy spectrum at t = 17.4 for £fc max = 1.48, 2.97 and 6.01; 
the dashed black line indicates k 2 power-law scaling, c) Tem- 
poral evolution of kinetic energy spectrum; the solid red lines 



correspond to fits using Eq.(62l and the dashes black line 
indicate k~ 3 power-law scaling, d) Temporal evolution of ef- 
fective self-truncation wavenumber k, 
resolutions. 



(Eq.(62l) at different 



enough? Note that this problem is related to the classical 
Fermi-Pasta-Ulam-Tsingu problem [32]. 

To try to answer this question within the Taylor-Green 
framework would be computationally very expensive as 
long runs should be performed at arbitrarily high resolu- 
tions. 

A simple alternative idea to study this problem is to 
use initial data for the TGPE generated by the SGLE 
with a variable truncation wavenumber fc™, set to a tar- 
get value of k c , smaller than the maximum truncation 
wavenumber fc max allowed by the resolution. This SLGE- 
generated initial data can then be used to run the TGPE 
at a given value of £k c with arbitrarily large values of 
£&max- A number of runs were performed at resolution 
64 3 with various values of fc™, £, and initial energy e ln 
(see legend on Fig 7]) . The result of these computations 
are compared with the above Taylor-Green runs (see 
Figj6]) and displayed on Fig{7j Because of the steep de- 
cay of the energy spectrum for k 3> k Cl the self-truncation 



k 2 e^ n (k)dk 
\ 3 J fcmaX e ^n(k)dk 



5 Jo 
3 



(64) 



A general growth in time of k c is apparent on FigjTja 
for both the Taylor-Green runs and the SGLE-generated 
initial data, showing similar behavior. In order to check 
for self-similar regime a parametric Log-Log represen- 
tation dk c /dt v.s. k c has been used on Figj7jb and 
FigjTjc. With this representation, a self-similar evolution 
k c (t) ~ f corresponds to a line of slope \ = i 7 ! ~ I)/ 7 ?- 
Figure [7}b, shows transient self-similar evolutions, that 
all terminate by a vertical asymptote, corresponding to 
logarithmic growth (77 = 0). This self-truncation takes 
place for small values of k c /k mllK strongly suggesting that 
the self-truncation happens in a regime independent of 
cut-off. Finally, Fig[7]c suggests that, depending on ini- 
tial conditions, self-truncation can take place at arbitrar- 
ily values of £fc c . 

As the dynamics of modes at wave-numbers larger than 
k c is weakly nonlinear, it should be amenable to a descrip- 
tion in terms of wave turbulence theory; this could per- 
haps explain the slowdown of the thermalization in this 
zone. The new regime indicates that total thermaliza- 
tion is delayed when increasing the amount of dispersion 
(controlled by £fc max ) but is preceded by a partial ther- 
malization (quasi-equilibrium up to k c ) within a PDE. 

We now turn to estimations of order of magnitude rel- 
evant to physical BEC. At low-temperature, the GPE is 
known [3] to give an accurate description of the (classical) 
dynamics of physical BEC at scales larger than the in- 
teratomic separation I. At finite temperature the TGPE 
gives a good approximation of Bose-Einstein condensate 
(BEC) only for the phonon modes with high occupation 
number, see [31 US]- At very low temperature thus only 
a limited range of low-wavenumber density waves are in 
equipartition. 

This limited range has consequences on the low- 
temperature thermodynamics of BEC that can be ob- 
tained by the following considerations. The equiparti- 
tion range is determined by the relation k < fc cq with 
Hu)B(k cq ) — ksT where the Bogoluibov dispersion rela- 
tion wb(A;) is given by 



uj B {k) 



9\Ao\ 



Am 



-k 2 . 



(65) 



The coherence length £ defined in Eq.([7]) can be expressed 
in terms of the s-wave scattering length a defined by g = 
Ait ah 2 Jm and the mean inter-atomic particle distance t = 
n- 1 / 3 « |A,|- 2 / 3 as 



-1/2 



1 



8ir \a 



1/2 



(66) 



For weakly interacting BEC the coherence length thus 
satisfies £ 3> £■ The equipartition wavenumber fc eq ex- 
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FIG. 7: (Color online) a) Evolution of the self-truncation wavenumber k c . Curves i-iv: £ = 2V2/5, K n = 4, e in = 0.1, 0.2 , 0.4 , 1; 
v: £ = V2/10, kf = 8, e in = 0.2; vi-viii: £ = ^2/5, e in = 0.1, 0.2,0.4 (i-viii in resolution 64 3 ); ix-xi: Taylor-Green resolutions 
64 3 , 128 3 and 256 3 . b-c) Parametric representation dk c /dt v.s. fc c /fc max and dk c /dt v.s. £fc c (same labels as in Fig. a). 



plicitly reads 



k — 

^eq — 



v /32fc|m 2 T 2 C 4 + h 4 1 
4^ 4^2 



1/2 



(67) 



Using the Bose-Einstein condensation temperature of 
non-interacting particles (valid for a <C £) [3] 

i 2/3 



2irh 2 



C(l) 



(68) 



where C(3/2) = 2.6124 . . ., the equipartition wavenumber 
k eq can be expressed as 



1 

2£ 



i28tt 2 £ 4 t 2 



1 1/2 



(69) 



Observe that fc eq varies from fc eq = at T = to 
wavenumber of order fc oq ~ i~ l at T\ and it is equal 
to fc 4 = 2tt/£ at T* defined by 

T* =47rv / l + 8^ 2 C(^ 2/3 ^T A . (70) 

Thus the thermodynamic of physical BEC at low- 
temperature (e.g. specific heat scaling as T 3 ) can be 
recovered from the TGPE thermodynamics by setting 

In experimental turbulent weakly interacting BEC 
such as [13] the value of £fc oq is large because T*/T\ ~ 
I <C 1 and therefore the corresponding TGPE should 
have a large £/c max . Thus the thermalization slowdown 
caused by the dispersive bottleneck should be in princi- 
ple present in physical BEC, unless it is overwhelmed by 
other relaxation mechanisms |4"4"1. 



V. METASTABILITY OF COUNTERFLOW, 
MUTUAL FRICTION AND KELVIN WAVES 

Counter-flow states with non-zero values of momen- 
tum generated by the new SGLE algorithm and their 



interaction with vortices are investigated in this section. 
The counter-flow states are shown to be metastable under 
SGLE evolution; the spontaneous nucleation of vortex 
ring and the corresponding Arrhenius law are character- 
ized in section |V A| Dynamical counter-flow effects are 
investigated in section |VB| using vortex rings and vor- 
tex lines patterns that are exact solutions of the GPE. 
Longitudinal and transverse mutual friction effects are 
produced and measured. An anomalous translational ve- 
locity of vortex ring is exhibited and is quantitatively 
related to the effect of thermally excited finite-amplitude 
Kelvin waves. Orders of magnitude are estimated for the 
corresponding effects in BEC and superfluid 4 He. 



Metastability of grand canonical states with 
counterflow 



1. Thermodynamic limit of states with nonzero counterflow 

The counterflow states with W ^ are determined by 
thermal fluctuations around the minima of the energy F 
Eq. ( 23 ) . These minima correspond to the solution of 



SF 







2m 



V 2 ip+gV G [\ip\ 2 ]ip-nip+ihW-Vip (71) 



that are plane-waves of the form 



[i — mW • v s 



(72) 



where the velocity v s indexes the different solutions. 

In the thermodynamic limit, the Galilean group de- 
fined by the transformations ( 12p5 1 is continuously in- 
dexed by the velocity vq. All wavefunctions (72) are 



thus equivalent by Galilean transformation (and redef- 
inition of the chemical potential). Under the Galilean 
transformation ( 12 ) the energy F is transformed as F' — 
F — (raW • vq — mv G /2)N + vq • P. Note that, among 
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all the minima of F the one with v s = W minimizes F' . 
This state corresponds to a condensate moving with uni- 
form velocity W. The W • P term is thus only imposing 
a Galilean transformation of the global minimum. 

However, when working in a finite volume, the Galilean 
transformation is quantized (see Eq.(19)). The minima 
of F' of lowest energy then corresponds to a condensate 
moving with the quantized uniform velocity v s that is the 
closest to W. At finite temperature and volume, when 
W is not too large with respect to the velocity quantum 
in Eq.(19l, we have two ways to produce momentum in 
the system. The first one corresponds to Galilean trans- 
formations: v s 7^ in (72). The second one to fluctu- 



ations of the exited phonons, with v s = in ( 72 ) and 
the momentum of phonons imposed by the term W ■ P 



in the grand canonical distribution ( 22 1 . Metastability is 
thus expected when W ^ with quasi-equilibrium corre- 
sponding to condensates at different wavenumbers with 
an energy barrier between each of those states. 

In the context of the Landau two-fluid model Q] the 
velocity v s of the condensate corresponds to the super- 
fluid velocity and the momentum carried by the exited 
phonons is written as P = p n ( v n — v s) where p n and v n 
are called the normal density and velocity respectively. 
The counterflow velocity defined by W = v n — v s is a 
Galilean invariant. 

The above discussion shows that, in general the vari- 
able W in the SGLE (42) corresponds to W = v n . In 
the thermodynamic (infinite volume) limit W = v s and 
there is thus no counterflow W = v n v s = 0. For 
finite-size systems, in general v s ^ W and W ^ 0. 

We thus define (when v s — 0) the normal density by 
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FIG. 8: (Color online) Counterflow dependence of momentum 
Pz (w x = W y =0). Inset: Temperature dependence of p n = 
|^ . b) Histograms of momentum P z and — P z (in log — 

lin) with no counterflow at T = 1. No asymmetry is observed, 
c) Histograms of momentum P z and —P z with counterflow 
W z = .4 at T = 1. An asymmetry, induced by counterflow, 
is apparent. Observe that both histograms are centered at 
P. = 0. 



but the non-zero counterflow induces an asymmetry in 
the statistical distribution that yields a non-zero value 
for the mean momentum. 



w z =0 



Thermodynamics of metastable states at small 
temperature and small counterflow 



3. Spontaneous nucleation of vortex rings and Arrhenius 

law 



In order to first validate the SGLE in the presence of 
counterflow two scans are performed at constant density 
using a resolution of 64 3 and £fc max = 1.48. The con- 
densate is set at k = in the SGLE initial data and 
the temperature is fixed to T = 0.2. This low temper- 
ature allows us to increase the value of the counterflow 
w z (hereafter we set w x = w v — 0) keeping the con- 
densate at k = 0. The dependence of the momentum 
P z on w z is presented in FigMa. The solid line corre- 



sponds to the low-temperature calculations (Eq. ( 59 ) and 
appendix [b] Eqs.(B9)). The second run correspond to 
a temperature scan (at low counterflow w z = .1). The 
temperature dependence of p n is displayed together with 
the low-temperature calculation on the inset of Fig|S]a. 

Figure[8]a-b display histograms of P z and — P z in phys- 
ical space, both obtained at T = 1 with the condensate 
at k = but with zero and non-zero counterflow. Ob- 
serve that the histograms are both centered at P z = 



At temperatures and counterflow velocities large 
enough the stochastic process defined by the SGLE can 
jump between the different metastable states discussed 
above in section |V A 1| In this section, we show how 
the different states are explored, under SGLE evolution, 
by spontaneous nucleation of vortex rings. To wit, we 
present a numerical integration of SGLE at resolution 
64 3 with £/c max = 1.48. With this choice of parameters 



the velocity quantum (19) is fixed to 0.2. The tempera- 
ture is set to T — 0.775 and the counterflow to w z = 0.8. 
The condensate is set at k = in the SGLE initial data 
and the density is kept constant to p = 1. 

The temporal evolution of the momentum P z is dis- 
played in Figj9]a (right scale). Observe that the system 
first spends some time at the state (/) with P z ss 0.05 
and that, around t = 55, it jumps to the state (II) with 
P z w 0.225. These two metastable states correspond to 
quasi-equilibrium at k = and k = 1 as is apparent in 
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FIG. 9: (Color online) a) Temporal evolution of |Aoj 2 and 
\Ai 2 (left scale) under SGL dynamics. Observe that there 
are two quasi-stationary states (7 and 77) and the conden- 
sate makes a transition from k = to k = 1. The temporal 
evolution of the momentum p z is displayed in the same plot 
(right scale). Observe that transition from one state to the 
other is accompanied by an increase of momentum, b) 3D 
visualization of density at t — 54.5, t — 56 and t = 60.5; the 
grey (blue) clouds corresponds to density fluctuations and the 
vortices are displayed as grey (red) isosurface (see colorbar on 
Fig 12 below), c) Histogram of momentum p z and —p z at the 
two quasi-stationary states (7 and 71) in lin — lin plot, d) 
Arrhenius law : data form SGL dynamics (points) and theo- 
retical Eq.{74| (solid line). 



Fig|9]a (left scale) where the temporal evolution of |^4n| 2 
ariclf-Ail 2 (see Eq.(16)) are displayed. 

In order to illustrate the dynamic of the condensate 
jump from k = to k = 1 we now present 3D visualiza- 
tion of the density at t = 54.5, t = 56 and t = 60.5 on 
Fig|9]b. To produce this figure, the wave-function if) was 
first low-pass filtered and the density was then visualized 
using the VAPOR software. At early times (t < 50, pic- 
tures not shown) no vortices are present in the box. At 
t ~ 54 a vortex ring is nucleated. It then increases its 
size under SGLE evolution until it reconnects with the 
neighbor rings (recall that periodic boundary condition 
are used). The ring finally contracts and disappears (pic- 
tures not shown). During this evolution, the local phase 
defect of the ring becomes global and changes the con- 
densate wavenumber. Histograms of momentum P z and 
— P z in the two metastable states 7 and 77 are presented 
on Fig|9]c. Observe that both momentum histograms of 
metastable states are asymmetrical (as it was the case 



on Figjljc). However, note that I is centered at P z — 
and 77 at P z — 0.2, respectively corresponding to the 
wavenumbers k = and k = 1. 

It is well known that the escape time of a metastable 
quasiequilibirum is given, in general, by an Arrhenius law 
031461 



to 



-f3AF 



(74) 



where AF is the activation energy of the nucleation so- 
lution and t c is a characteristic time. Here, the nucle- 
ation solution is given by a vortex ring that satisfies 
= 0. The energy barrier is thus determined by 
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expressions for the energy 77 r 
and the radius are given by 
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where is the density at the infinity and a is a core 
model-depending constant with value a = 0.615 for the 
GPE vortices [4]. Formulae (75 77) and the value of a 
have been numerically validated in reference [47] using a 
Newton method 



In order to numerically check that the escape time in- 
deed follows an Arrhenius law we now perform runs with 
with £fc max = 1.48 and resolution 32 3 . The counter-flow 
is fixed at w = 1.4 and the condensate is set initially 
at k = (constant density p = 1). At each fixed tem- 
perature T, several numerical integration of SGLE are 
performed and the escape times for the condensate to 
leave the wavenumber k = are measured. These es- 
cape times are then averaged over more than 10 realiza- 
tions. Figure [9jd displays the escape time t csc obtained 
in this way as a function of the inverse temperature 1/T 
in log — lin. The slope of the solid line is computed using 
the analytic formulae ( 75|77 ) of AF. Both, numerical 
and theoretical Arrhenius laws are in good agreement. 
The main consequence of this Arrhenius law is that it 
is practically possible to use the SGLE dynamics to pre- 
pare metastable states with finite value of counterflow 
and lifetime quantitatively given by (74). 



B. 



Dynamical effects of finite temperature and 
counterflow on vortices 



We now turn to the study the dynamical effects of 
counterflow on TGPE vortex evolution. To wit, we set 
up finite temperature and finite counterflow initial states 
that also contain vortices. Two cases are investigated: (i) 
vortex lines, in a crystal-like pattern that does not pro- 
duce self induced velocity and (ii) vortex rings, producing 
self induced velocity. 
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1. Lattice of vortex lines 

To numerically study the effect of counterflow on vor- 
tices we prepare an initial condition lattice consisting in 
a periodical array (of alternate sign) straight vortices. 
This initial condition is the 3D extension of that used in 
ref.[5T] to study the scattering of first sound in 2D. The 
lattice is obtained using a Newton method [48-50J. It is 
an exact stationary solution of the (periodic) GPE. As 
the vortices are separated by a fixed distance d = it, they 
can be considered isolated in the limit £ <C d. Let us re- 
mark that this limit is automatically obtained when the 
resolution is increased at constant £fc max . To include tem- 
perature effect we prepare absolute equilibria ip cq using 
SGLE with the counterflow aligned with an axis perpen- 
dicular to the vortices in lattice- The initial condition 
i> = ^lattice x "0oq is then evolved, using the TGPE. The 
counterflow induces a motion of the lattice as is apparent 
on the 3D visualizations of the the time evolution of the 



component of the vortex filament to the counterflow for 
T = 0.5, 1 and w z — 0.4. The trajectories are obtained 



density that arc displayed on Fig 10 (Figure obtained in 
the same way as Figj9)b). Several runs were performed 




FIG. 10: (Color online) 3D visualization of density at t = 0, 
40, 60 and 120 at temperature T = 1 and counterflow W = 
0.4. The grey (blue) clouds correspond to density fluctuations 
and the crystal-like vortex lattice is displayed in grey (red) 
isosurfaces 



at different resolutions (with £fc max = 1.48 ), t emperature 
and counterflow values (see legend on Fig[TT|b). 

Figure [ll]a displays the temporal evolution of 
(R\\,R±) the respectively parallel and perpendicular 




FIG. 11: (Color online) a) Trajectory of a straight vortex 
in the crystal pattern for T = 1, T = 0.5 and w z = .4, at 
resolution 64 3 . Inset: run with T = 1 until t = 600. b) Tem- 
perature dependence (T\ = 3.31) of the advection velocity 
v«/w z for the lattice and Adl/ui for the vortex rings (reso- 



lutions 32 3 -128 3 ). The dashed line corresponds to to Eq. ([79 1 
with B' = 0.83 and the solid line to the theoretical prediction 
(901. c) Temporal evolution of the square of the length of 
the vortex ring for different values of counterflow, T = 1 and 
initial radius R = 15£ at resolution 64 3 . 



by first averaging along the direction of the vortices, then 
the (averaged) coordinate of the vortices is found by seek- 
ing the zero of the reduced 2d wavefunction. Observe 
that the vortex, originally located at (^, 4^), moves in 
the direction of the counterflow and its velocity clearly 
depends on the temperature. It is apparent that a per- 
pendicular motion is also induced at short times. This 
motion has two phases, the first one is related to an adap- 
tation and makes the crystal-like lattice slightly imper- 
fect. Then the perpendicular motion almost stops (a very 
small slope can be observed for long time integration). 
The initial phase where the parallel and perpendicular 
motions have similar velocities lasts longer when £/d is 
decreased by increasing the resolution (data not shown) . 
Observe that the imperfection of the lattice in the final 
configurations is almost equal for the two temperatures 



presented in Fig 11a, but the parallel velocities are con- 
siderably different. Thus, the self-induced parallel veloc- 
ity caused by the slight imperfection of the lattice is very 
small and is not driving the longitudinal motion. 

We now concentrate on the measurement of Rii for 
which the present configuration is best suited. 

Ru has a linear behavior, that allows to directly mea- 
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sure the parallel velocity V\\. The temperature depen- 
dence of vu/w z is presented on Fig 11 b for different val- 



ues of w z and c£/£ (corresponding toThc different resolu- 
tions) . 

For superfluid vortices the standard phcnomcnological 
dynamic equation of the vortex line velocity vl is [I] 



vl 



v sl + as' x (v n -v s) )-aV x [s' x (v n - v sl )], (78) 



where s' is the tangent to the vortex line, v s i is the lo- 
cal superfluid velocity (the sum of the ambient superfluid 
velocity v s and the self-induced vortex velocity Ui) and 
v n = w + v s is the normal velocity. The constants a, a' 
depend on the temperature. Let us remark that the exis- 
tence of the transverse force (related to the third term of 
r.h.s. in Eq 78 ) has been subject of a large debate in the 



low-temperature community in the last part of the 90's 
[52Tf58] and this controversy is still not resolved. Applied 
to the present case, Eq.(78) predicts v± — —aw z and 
vii — a'w z . The value of the constant a', related to the 
transverse force, depends on the normal density and the 
phonon-vortex scattering section. It can be expressed as 



B 



2p 



where B' is an order one constant [J|. A fit to the mea- 
sured values of wj|/io 2 yields B' — 0.8334, see FigfTTjp. 
We thus conclude that finite-temperature TGPE coun- 
tcrflow effects measured on i?y for the crystal pattern are 
in quantitative agreement with standard phenomenology 
(Eq.(78l). We have seen above that the effect on R± is 
of the same order of magnitude that the one on i?y but 
only in the initial phase as long as crystal imperfection 
does not come into play. 



finite-difference scheme version of the TGPE that exactly 
conserves the energy and particle number [30) . 

The temporal evolution of the square of the vortex 
length of a ring of initial radius R = 15£ at temper- 
ature T = 1 and counterflow w z = 0, 0.2 and 0.4 is 
displayed on Fig 11 c. For w = 0, the dynamics under 
TGPE evolution reproduces the Berloff ring contraction 
[5U] . The temperature dependence of the contraction ob- 
tained for w = (data not shown) quantitatively agrees 
with Berloff and Youd's results. A dilatation of vortex 
rings is apparent on Fig |ll| c for w larger than the mea- 
sured vortex ring velocity vj, = 0.23. 

However, vt, has a very strong dependence on temper- 
ature that is also present for w = 0. The temperature 
dependence of Aw^/u; where Avr, — — vl and is dis- 
played on FigjTTJb. We have checked that the velocity 
vl directly measured at T = is indeed given by u;. 
Equation |78| ) predicts (in the absence of counterflow) a 
translational velocity for the vortex ring vi, = (1 — a')ui. 
Observe that A^/it; is one order of magnitude above 
the transverse mutual friction coefficient measured on the 
crystal-like lattice. Note the presence of a large spread of 
the low temperature data for the ring (see the leftmost 
datapoints on Fig 11 b corresponding to T/T x = 0.04). 
At very low temperatures the effect is very weak. Thus, 
the corresponding measured values of Avr,/ui are influ- 
enced by errors on the measurement of position and ve- 
locity of the vortices that are caused by the finite-size 
of the mesh. In the future, these low temperature un- 
certainties on the determination of Avl/ui could be re- 
duced by performing runs at higher resolution together 
with more accurate (sub-grid) measurement of the vortex 
position. 



2. Vortex rings 

We now turn to study the effect of counterflow on vor- 
tex rings. The initial condition is prepared as in the 
previous section but with the lattice f/'iatticc replaced by 
a vortex ring V'ringj that is an exact stationary (in a co- 
moving frame) solution of GPE. The plane containing the 
vortex rings of radius R is perpendicular to the counter- 
flow and the rings are numerically obtained by a Newton 
method gHHSD]. 



In the case of vortex rings the general formula (78) 
yields 

R = -a(ui - w z ) (80) 

vl = v s + (1 - a')ui + a'w z , (81) 

where u; denotes the ring velocity at zero temperature, 
explicitly given by Vri n g in formula ( 75 1 (replacing R* 
by the corresponding radius). In the special simple case 
w z — 0, a finite-temperature contraction of the vortex 
ring is predicted. This transverse effect effect was first 
obtained and measured by Berloff and Youd, using a 



3. Anomalous translational velocity and Kelvin waves 

In this section we relate the finite temperature slow- 
down (see the top line of Fig 
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b) to the anomalous 
translational velocity of vortex ring with finite-amplitude 
Kelvin waves that was reported in refs. (SnUfSJj. Indeed, 
Kelvin waves are clearly observed in 3D visualizations of 
vortex rings driven at finite-temperature by the TGPE 
as it is apparent on Fig |12| (obtained in the same way 
that Fig|9]b). 

Following reference |60| . Kelvin waves of amplitude 
A and wavelength 2-kR/N on a ring of radius R are 
parametrized, in cylindrical coordinates r, <j) and z, as 



x = (R + A cos N(f) cos 4> 
y = (R + A cos Ncf>) sin </> 
z = — As'md). 



(82) 
(83) 
(84) 



In the limit N 3> 1 the dispersion relation uj(k) of the 
Kelvin wave ( 82"|84 ) is given by [59, 
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FIG. 12: (Color online) 3D visualization of density at t — 18, 
19, 20 and 21 at temperature T = 1. The grey (blue) clouds 
correspond to density fluctuations and a vortex ring of radius 
R = 20£ with thermally excited Kelvin waves is displayed in 
grey (red) isosurfaces 



which yields the value of A 2 N 2 / R 2 as function of T: 
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ir 2 p c 



m 2 kBT 
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The equipartition law 
as the classical limit of t 
puted by Barcghi et al. 



can also be directly obtained 
le quantum distribution com- 
61], up to a redefinition of the 
core constant model a (see Eq.(25) in reference [ST]). Let 
us remark at this point that, at low temperature and in 
non-equilibrium conditions, the presence of a Kelvin wave 
cascade at scales between the inter-vortex distance and £ 
(see references 62,63 ) can lead to a different dependence 
of the amplitude on the wavenumber. 

We finally assume that the slowing down effect of each 
individual Kelvin wave is additive and that the waves 
populate all the possible modes. Kelvin waves are bend- 
ing oscillations of the the quantized vortex lines, with 
wavenumber k < 27r/£. The total number of modes can 
thus be estimated as 



Mo 



lvin 



2irR/£. 



(89) 



Replacing A 2 N 2 /R 2 in Eq.(|86j by Eq.(|88j) and multi- 
plying by the total number of waves A/jceMn we obtain 
the following expression for the anomalous translational 
effect due to thermally exited Kelvin waves 



2k B Tm 2 



1 



Ui 



TTPoo^ 2 log m 



(90) 



where k = N/R and a is the core model-depending con- 
stant in formula (75). 



The anomalous translational velocity caused by an ex- 
cited Kelvin wave was first reported by Kiknadze and 
Mamaladze [SS] in the framework of the local induction 
approximation (LIA). The effect was then obtained and 
numerically characterized within the Biot-Savart equa- 
tion by Barenghi et al. [60] . The anomalous translational 
velocity u a of a vortex ring reads (in the limit N ^> 1, 
see Eq.(26) of reference [59"] ) 



, A 2 N 2 . 



(86) 



where it; = Vring is the self-induced velocity ( 75 ) without 
Kelvin waves. 

The variation of the energy of a vortex ring caused by 
a (small amplitude) Kelvin wave can be estimated as 



AE _ Coring AL 

dR 2ir 



(87) 



where H T i ng is the energy given by Eq.(77) and the 



length variation AL produced by the Kelvin wave (82 
84 1 is given, at lowest order in the amplitude A/R, by 
AL = ttA 2 N 2 / R. Assuming equipartition of the energy 
of Kelvin waves with the heat bath implies AE = ksT, 



The temperature dependence of the equipartition esti- 



mate (90) of the thermal slowdown is plotted on Fig 11 b 



(top straight line) . The data obtained form the measure- 
ments of the rings velocity in the TGPE runs is in very 
good agreement with the estimate ( 90 ) . 

As discussed in refs. [3] [T3] the TGPE gives a good 
approximation to physical (quantum) Bose-Einstein con- 
densate (BEC) only for the modes with high phonon oc- 
cupation number. In this sprit quantum effects on the 
Kelvin waves oscillations must also be taken into account 
to obtain the total slowing down effect in a BEC. The 
TGPE estimation (90) can be adapted to weakly inter- 



acting BEC by the following considerations. 

At very low temperature, because of quantum effects, 
only a limited range of low- wavenumber Kelvin waves are 
in equipartition. This range is determined by the relation 
k < fc cq with huj(k cq ) = k^T and the dispersion relation 



(85), it reads: 



kftT 2m 

8R\ 



- y^[ln( ? 
and can also be expressed as 



) 



(91) 
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/ 4-7T n 2 / 3 
C(§)V3[ln(ip )- fl ] 



1/2 



(92) 
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where T\ is the Bose-Einstein condensation temperature 



of non-interacting particles ( 68 ) and the relation between 
the interatomic distance £ and the vortex-core size £ are 



given in Eq. ( 66 ) 

Observe that fc, 
wavenumber of order k, 



cq varies from fc oq = at T = to 
t~ x at T\ and it is equal to 



eq 



k£ = 2tt/£ at T* defined by 



T*=8vr 2 C(f) 2 / 3 [My) 



a . 

£ 



(93) 



Therefore at temperatures T* < T <T\ the energy of all 
Kelvin waves are in equipartition and equation ( 90 ) thus 
applies directly. 

It is natural to suggest that an additional effect, caused 
by the quantum fluctuations of the amplitudes of Kelvin 
waves, will take place at low temperatures T < T* . This 
quantum effect can be estimated by using the standard 
relation for the energy of the fundamental level of a har- 
monic oscillator AE — fuo(k)/2. Applied to the Kelvin 
waves, this relation yields the fc-independent quantum 
amplitude Aq = m/4ir 2 Rp. The quantum effect can thus 
be estimated as the sum 



Mum, A 2 qN 2 
^ B 2 



clvin 



2m?r 647T 5 / 2 



3i? 2 



3^2 



- X 3/2 



(94) 

The total effect is obtained superposing the thermal effect 
and the quantum effect and the final result is 
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Ui 
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(95) 



(96) 



where C[R/£] = log ~ a. 

In the case of superfluid Helium, where a ~ £, the 
GPE description is only expected to give qualitative pre- 
dictions and, at best, order of magnitude estimates (see 
ref.[4]). It is thus difficult to extend the above consider- 
ations, obtained in the case of weakly interacting BEC 
with a <C £, to Helium. 

Nevertheless the results obtained above in the weakly 
interacting case strongly suggest the presence of new 
slowing down effects, not included in the usual mu- 
tual friction descriptions of Helium that predicts ~ 
Pn/p ~ {T/TxY. The new effects, because of their tem- 



perature dependence (see Eq.(96)), should be dominant 
at low-temperature. 

The zero-temperature quantum slowdown is indepen- 
dent of the ring diameter and the finite temperature ef- 
fects are stronger for small rings. Time of flight measure- 
ments of vortex rings in 4 He could be used to determine 
the translational velocity. The effect could also be stud- 
ied in ultra-cold atomic gases BEC. For these systems 
the effect of the inhomogeneity of the superfluid should 
be taken into account 



VI. CONCLUSIONS 

In summary, our main results were obtained by making 
use of a stochastically forced Ginzburg-Landau equation 
(SGLE) that permits to efficiently obtain and control 
truncated Gross-Pitaevskii absolute equilibrium. This 
allowed us to show that the condensation transition ob- 
served in references [13j[27l[28] corresponds to a standard 
second-order transition described by the A</> 4 theory. 

We also found that thermodynamic equilibrium can be 
obtained by a direct energy cascade, in a way similar to 
that of Cichowlas et al. [2T] , accompanied by vortex anni- 
hilation as a prelude to final thermalization. Increasing 
the amount of dispersion of the system a slowdown of 
the energy transfer was produced inducing a partial ther- 
malization independently of the truncation wavenumber. 
This new thermalization regime opens up an avenue to 
a further investigation of vortex dynamic in co-flowing 
finite-temperature superfluid turbulence. In this con- 
text it would be interesting to study in the future, us- 
ing a much higher resolution than in the present work, 
the dispersive bottleneck. In particular to investigate 
the possibility of the coexistence of a well established 
turbulent Kolmogorov cascade followed by a dispersive- 
induced partial-thermalization zone. 

Using the SGLE in the presence of a counterflow we 
observed that the counterflow can block the contraction 
of vortex rings reported by Berloff and Youd |30j and 
also induce a dilatation. We directly measured the mu- 
tual friction coefficient related to the transverse force. 
An unexpected result was found by immersing a vortex 
ring in a finite-temperature bath: a strong dependence 
of the translational velocity in the temperature was ob- 
served. This effect was an order of magnitude above 
the transverse mutual friction effect. We explained this 
effect by relating it to to the anomalous translational 
velocity due to finite amplitude Kelvin waves that was 
previously found by Kiknadze and Mamaladze [SS] and 
Barenghi et al 60J. Assuming equipartition of the en- 
ergy of the Kelvin waves with the heat bath yields a 
formula that gives a very good quantitative estimate of 
the numerically observed effect. This new formula also 
gives an experimentally-testable quantitative prediction 
for the thermal slowdown of vortex rings in weakly in- 
teracting Bose-Einstein condensates and superfluid 4 He. 
In this context, it would be interesting in the future to 
study (using a higher resolution than in the present work) 
the vortex dynamic of counter-flowing finite-temperature 
superfluid turbulence. Note that, in the context of BEC 
(where experiments are performed within a confining po- 
tential) the wavefunction tp must be expanded using an- 
other basis of orthogonal independent functions than the 
Fourier modes (e.g. the eigenfunctions of the harmonic 
oscilator) , see ref. |65j where the apparent arbitrariness of 
the truncation parameter is also discussed. About this 
last point, also see the discussion around Eqs.(65) and 
( 70 1 at the end of Sec. IV B of the present work about 



the upper limit fc oq of the equipartition range that follows 
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from the quantization of phonons in a physical BEC. 

The TGPE dynamics was thus found to contain many 
physically sound phenomena of finite-temperature super- 
flows. This strongly suggests the possibility to obtain 
the propagation of second sound waves in the TGPE. 
Some preliminary results support this conjecture (data 
not shown), however very high resolutions seem to be 
needed and this will be the subject of a future work. 



Acknowledgments 

We acknowledge useful scientific discussions with G. 
During and S. Rica. The computations were carried out 
at IDRIS (CNRS). 



Appendix A: Conservation Laws and Dealiasing 

In the standard incompressible Euler case, for 
quadratic nonlinearities and quadratic invariants, the 
system can be correctly dealiased using the 2/3— rule that 
consists in truncation for wavenumber |k| < fc max = N/3, 
where N/2 is the largest wavenumber of the discrete sys- 
tem. With this procedure, one third of the available 
modes are not used. Such discrete dealiased pseudo- 
spectral system exactly conserve the quadratic invariant 
and is therefore identical to the original Galerkin trun- 
cated system. 

In the TGPE case, the problem is more complicated 
because the equation is cubic and the invariants are quar- 
tic. Let us first recall Parseval's theorem that states 
/ d 3 xf(x)g*(yi) = VJ2kfk9t, where / k and g k are the 
Fourier transform of / and g. This identity remains valid 
in truncated systems and it holds whether the functions 
are dealiased or not. The integration by parts formula is 
a consequence of Parseval's theorem: 



3 k 



d 3 x 77— 3 



Remark that the product rule (fg)' = f'g + fg' is only 
valid if the fields are dealiased. 

The conservation of the total number of particles is 
directly obtained using the GPE 



dN 
~dt 



ih 
2m 



where the last equality is a consequence of the Parse- 
val identity and is thus true independently of dealiasing. 
Similar relations lead to the conservation of the energy 
H. 



Using the dealiased TGPE ( 18 1 the conservation law 
for the momentum reads 



dt 



= 2g / d 3 x [(^ G [M 2 ]) H 2 +P G [|Vf]^#| 2 



If ip is dealiased the 2/3— rule implies that 



d 3 x{V G [\^\ 2 ]^)d^ = / d 3 xV G [V G M 2 }4>] drf 



d s (v G [\M 2 ]i>) = {d 3 v G [\M 2 ])^ + v G U\ 2 ]d^ 

dP 

it follows that —rf- = 0. Without a Galerkin projector 



in Eq.(18l the aliased field would obey {\ip\ 2 ijj) djtp + 
(\ip\ 2 ip) dj'tp dj(\ip\ 4 ) and the conservation of momen- 
tum would therefore be lost. 

Conservation of N, H and P can be numerically 
checked by using absolute equilibria with non-zero mo- 
mentum. The conservation of P is ensured only if the 
system is dealiased. The error of aliased runs grow up 
to a 50% in a few units of time and is independent of 
the time-step (data not shown). We thus believe that it 
would important to explicitly check the conservation of 
momentum when using finite-difference schemes, even if 
they exactly conserve the energy and the particle num- 
ber. 



Appendix B: Low-temperature calculation of 
thermodynamic functions 

We are interested in computing the grand partition 
function Z in Eq.(47) where F = H — fiN — W • P is 



written in terms of 



Fourier amplitudes as 



H 

V 

N = 
P, = 



E 



h 2 k 2 

~2m 



k 



k 3 +ki ^k 2 ^k 4 +k 2 ^k 3 ,-k 4 



(Bl) 
(B2) 



(B3) 



where — if k > fc max and the second sum in H is 
over ki, k 2 , k 3 , k 4 . 

The saddle-point is determined by the condition J^pr — 
fJ-oAoVSk o = which, separately written for k = and 
k ^ 0, explicitly reads 

(g\ A 1 2 -fi + ^)A + 2g £ A \A kl \ 2 (B4) 

ki#0 

+g ^ A *i A k-*i A -^ = 

ki,k 2 #0 

h 2 k 2 

A k - fiA k - HW ■ k A k (B5) 



2m 



'9 E ^ki^ka+k^k+ka = 



from which Eq. ( 48 ) follows 



(Al) 



To diagonalize F = H — fj,N — W • P we first apply the 
Bogoliubov transformation to H—fiN and then show that 
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P is also diagonal in this basis. Replacing B p , defined by 
the transformation (52 1, in H — fiN (recall that p = Kk) 



and then imposing the diagonalization determines the 
coefficient L TI 



-*p- 



Lp — 



-2|A |V 



H + e(p) 



\Ao\ 2 9 



where e(p) is given by 



e(p) 



2|A | 2 i 



2m 



/' 



\A \ 4 9 2 - (B7) 



The dependence of the entropy on the phase-space nor- 
malization constant is manifested by the presence of the 
logarithm term in S and A^- Notice that the function 
S + flXjsf is, however, completely defined. Also note 
that the pressure p must be computed, by definition, at 
constant total number of modes Af. All the thermody- 

(B6) 

namic relations discussed in section [Tl B| can be explicitly 
checked on the low-temperature expressions. The previ- 
ous formulae are represented as the solid lines that are 
confronted with the SGLE numerically generated data in 
Fig(2]a-b. 



The dispersion relation ( 54 1 is obtained by replacing 
|A 1 2 by its saddle-point value Eq.(48). 



We now express P in the Bogoliubov base. Using ( 52 ) 
directly yields 



lnJ 2 |BJ 2 - 



|B_ P | 2 + (u* p v* p B p B-p + c.c). 

(B8) 



Replacing Eq. ( B8 ) in the definition of P ( B3 1 , the last 



two terms vanish by symmetry and using the relation 



1, the momentum (B3) reads P 



£ P Pl-E?p| 2 ^- Formula (53) is then finally obtained by 
gathering H — fiN and W • P. 

The mean value of the condensate amplitude is ob- 

All the thermodynamic 



tained as V\A \ 2 = - 9n 



Mo=0 



variables are directly generated by first putting /io = in 
(p>9|) and then by differentiation, using relation (32). The 



fluctuations of the number of particles are computed as 
(57V 2 = — /3~ 1 f^- These quantities are explicitly listed 
below. 
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